Coupled-Cluster Approach to Electron Correlations in 
the Two-Dimensional Hubbard Model 
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We have studied electron correlations in the doped two-dimensional (2D) Hubbard model by using 
the coupled-cluster method (CCM) to investigate whether or not the method can be applied to 
correct the independent particle approximations actually used in ab-initio band calculations. The 
double excitation version of the CCM, implemented using the approximate coupled pair (ACP) 
method, account for most of the correlation energies of the 2D Hubbard model in the weak {U /t ~ 
1) and the intermediate U/t regions {U/t ~ 4). The error is always less than 1% there. The 
ACP approximation gets less accurate for large U/t {U/t ~ 8) and/or near half-filling. Further 
incorporation of electron correlation effects is necessary in this region. The accuracy does not depend 
on the system size and the gap between the lowest unoccupied level and the highest occupied level 
due to the finite size effect. Hence, the CCM may be favorably applied to ab-initio band calculations 
on metals as well as semiconductors and insulators. 



PACS numbers: 31.15.Dv,71.10.Fd,71.15.-m. 

Electron correlation in solids is the origin of various 
macroscopic quantum phenomena such as magnetism 
and high- Tc superconductivity. Interests in these 

fields are gradually shifting from materials which have 
simple chemical compositions and structures to more 
complicated ones. Besides the importance of simplified 
models such as the Hubbard and Heisenberg models in 
statistical physics, the importance of approximate many- 
body theories which can be applied to ab-initio calcula- 
tions is rapidly increasing. 

The coupled cluster method (CCM) has been recog- 
nized as one of the most successful approximate methods 
for quantum many-body problems in nuclear physics |^] 
and ab-initio quantum chemistry. |^ The method has also 
been successfully applied to the electron gas problem, |^] 
but there have been few applications to quantum lattice 
problems. Within the limited number of examples in the 
literature, we find that the method was successfully ap- 
plied to quantum spin models, Q the one-dimensional 
Hubbard and Pariser-Parr-Pople (PPP) model, Q and 
the two-dimensional (2D) Hubbard model at half-filling 
and with one hole. ||T^ These successful results suggest 
the CCM may be one of the most suitable methods to 
study electron correlations in solids as well as molecules 
and nuclei. 

To verify this expectation, we have made more rigor- 
ous tests of the CCM than the previous ones which were 
made on the one-dimensional Hubbard model on small 
clusters. Q We studied the doped 2D Hubbard model 
at closed shell fillings on finite size square lattices up to 
the 8x8 lattice by using the CCM. The total energies 
were compared with those obtained by using the pro- 
jector auxiliary field quantum monte carlo (PAFQMC) 
method and the the adaptive sampling quantum monte 
carlo (ASQMC) method. (An augmented version of 



the PAFQMC method which reduces the difficulty of the 
negative sign problem.) We have studied how the ac- 
curacy of the CCM depend on the Fermi degeneracy or 
magnitudes of the gaps of the system to be studied as 
well as magnitudes of Coulomb interactions and the elec- 
tron fillings. The former was not studied in the previous 
work, but it is very crucial if one tries to use the CCM 
for solid state electrons. 

We have studied the 2D Hubbard model: H = 
-tJ:i^,J)ai4aCJ. + H.C.) + Uj:^n,^n,^. We put i= 1. 

In t;he CCM, the exact totally symmetric nondegen- 
erate ground state wavefunction I^E*) is expressed on the 
basis of the exponential Ansatz: 

1*) =exp(T)|<i>), 

T — T„ , where r„ is a n-body excitation cluster 
operator. If $ is the Hartree-Fock (HF) wavefunction 
obtained in the Hartree-Fock limit (easily obtained for a 
lattice Fermion model), we may neglect the single excita- 
tion cluster Ti. The double excitation cluster may then 
be sufficient: T Ki T2 ■ This is an approximation called 
the coupled-cluster-double (CCD) approximation. Fur- 
ther incorporation of electron correlation is possible by 
taking into account of the remaining clusters r3,T4, • ■ •. 
It should be noted however that owing to the exponential 
form: exp(T), higher excited configurations than dou- 
bly excited configuration, such as quadruple excitation 
etc., are taken into account in the CCD approximation. 
If we use the normal product form of the Hamiltonian 
i?Ardefined in time-independent perturbation theory, jl^] 
the Schrodinger equation may be written as follows: 

{Hn - Ai;)exp(T)|$) 0, 

where /S.E = E — Ehf is the difference of the exact 
ground state energy and the Hartree-Fock energy. The 
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Schrodinger equation is solved by multiplying exp(— T) 
from the left hand side of the former equation and pro- 
jecting the resultant equation on excited configurations 
made from <!>: ^ 

($^/l[ffjvexp(T)]c|<i>) =0, 

where the subscript C indicates that only connected di- 
agrams are taken into accounted. The total energy E is 
given by: 

E = EhF + ($|[i/ArCXp(r)]c|$). 

The former equations are the coupled-cluster equations 
and they and the latter equation constitute a sufficient 
condition to the Schrodinger equation in the projected 
space. The advantage of the CCM is that exp(T) is ex- 
actly manipulated without artificial truncations owing to 
the presence of the subscript C, unlike other variational 
theories. This guarantees exactly the size consistency 
condition, which is difficult to satisfy in variational theo- 
ries. However, there is no variational upper bound given 
by the CCM energy. The double excitation cluster is 
expressed as: 

T2 = 1/4 E i'^^^^b'^J^ 

ijj.a.b 

where i and j are annihilation operators for single 
electron states occupied in $ and and are cre- 
ation operators for single electron states unoccupied in 
<i>. The CCD coupled cluster equation is given by: 
(e. + e, -ea- e^t^^ = (zj||a6) + 1/2 Ec<i(a^l 1^^)*?/ + 
y^EkiMklKl - EkdimjdKf - {bk\\zd)tl^ - 
{ak\\jd)tf, + {ak\\td)tf,) + EkUmcdMtt^f + Aitt^ - 
mqlt^^f + tp'g) - l/2{tpii + 4ft%f) + imftH] . The 
total correlation energy is given by: 

AE = l/4j2mab){t^^ + t^t';-t1t';). 

ijab 

{ij\\ab) is the anti-symmetrized two-electron integral de- 
fined by: {ij\\ab) = {ij\ab) — {ij\ba). is the i-th single 
particle energy. Each term in the CCD coupled-cluster 
equations can be assigned to a diagram, ft is known 
based on diagram theoretical arguments that the effect 
of the quadruple cluster T4 may be approximately taken 
into account by neglecting the 8, 9, 13, and 14-th terms in 
the right hand side of the CCD coupled-cluster equations. 
Such an approximation is known as the approximate cou- 
pled pair (ACP) approximation. |9|,[T^ 

We have calculated the total energy of the 2D Hub- 
bard model by using the ACP approximation. We have 
compared the ACP energies with exact energies obtained 
by using the PAFQMC and/or ASQMC methods. These 
are simply denoted by the QMC method, hereafter. The 
calculations were made on the 4 x 4, 6 x 6, and 8x8 lat- 
tices. Electron fillings p = Ne/Ng are such that electrons 
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FIG. 1. Total energies of the 4x4 Hubbard model cal- 
culated by using some approximations. Eg denotes the ex- 
act ground state energy calculated with the QMC method. 
p = 0.625. The energies are plotted as a function of U/t. 
Notations for the symbols are summarized in the legends. 

form closed shells on these finite size lattices, where Ne 
is the number of electrons and is the number of sites. 
We did not use the second order perturbative selection 
rules frequently used in ab-initio quantum chemistry cal- 
culations, which raise systematic errors in results when 
U/t is large. We have improved an existing fast algo- 
rithm by utilizing symmetries of two-electron integrals 
in the Hubbard model. It takes about 10 minutes of 
CPU time for a calculation of the 8x8 lattice on the 
Alpha workstation with the Alpha 21164/533 MHz CPU 
chip when we use the ACP approximation. It would not 
be difficult to do a CCM calculation even beyond a 128 
site lattice. We have not done this because the QMC 
calculations take much more CPU time than the ACP 
calculations. 

The ground state energies of the 2D Hubbard model 
on the 4x4 lattice calculated with various methods are 
plotted in Fig. |l]. The electron filling p is 0.625. The 
Restricted Hartree-Fock (RHF) method gives very poor 
values of the ground state energy even for very small 
values of U/t. The second order MoUer-Plesset (MP2) 
approximation works fairly well in the small U/t region 
{U/t ~ 1), but it fails in the intermediate {U/t ~ 4) and 
the large U/t regions {U/t ~ 8). The results obtained 
by using the ACP approximation are in good agreement 
with the QMC energies up to the intermediate U/t region 
on this lattice. The similar results obtained on the 8x8 
lattice are plotted in Fig. |[ p is 0.78125. Again, the 
RHF method is a very poor approximation even for very 
small value of U/t. The MP 2 approximation works only 
in the small U/t region. The ACP approximation works 
well up to the intermediate U/t region. The deviation of 
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FIG. 2. Total energies of the 8x8 Hubbard model calcu- 
lated by using some approximations. Eg denotes exact ground 
state energy calculated with the QMC method, p = 0.78125. 
The energies are plotted as a function U/t. Notations for the 
symbols are summarized in the legends. 

the energy calculated by using the ACP approximation 
from the energy calculated by using the QMC method on 
the 8x8 lattice is somewhat larger than that obtained 
on the 4x4 lattice when U/t = 8. 

The errors % of the ACP approximation defined by 
{~~Eacp + Eqmc)/Eqmc X 100 calculated on 4x4, 6x6, 
and 8x8 lattices are plotted against U/t in Fig. |[ p 
of the 4 X 4, 6 X 6, and 8 x 8 lattices are 0.6250, 0.7220, 
and 0.78125, respectively. In all cases studied, the errors 
are less than 1%, when U/t < 4. The errors grow as we 
increase U/t and they are larger than 2% when U/t = 8. 
The error in the 8x8 lattice reaches almost 7% when 
U/t = 8. Thus the ACP approximation breaks down in 
the large U /t region. The errors increase most rapidly 
in the case of the 8x8 lattice. The enhancement of 
the error in the larger lattice may be mostly due to the 
difference of p. p is closest to 1 in the 8x8 lattice in our 
cases. However, it is interesting to ask whether or not 
the enhancement is brought about by the decrease of the 
gap between the lowest unoccupied level (LUL) and the 
highest occupied level (HOL) in the finite size cluster. 

To clarify this point, we study the following gen- 
eralized Hubbard model: H — ~tJ2{i j)ai^lcr^jo- + 
H.C.) + UE.n.^n^i - E(,,,>(clc,. + H.C.) - 

*"E[ij]a(cLcjT + H.C), where (ij) and [i,j] denote 
the second and third nearest pairs of sites on the square 
lattice, respectively. We put t = l,t' — —0.2, and U = A. 
t" is a variable. We introduced anisotropy ±0.0001 on 
the X and y components of t and t" , as well as on the 
(1,1) and (1,-1) components of t'. We studied the 
8x8 lattice system with 78 electrons. If 0.0035 < <" 



FIG. 3. U/t dependencies of errors (%) of the ACP ap- 
proximation, calculated on 4 x 4, 6 x 6, and 8x8 lattices. The 
error (%) is defined by {-Eacp + Eqmc)/Eqmc x 100. p 
for the 4 X 4, 6 X 6, and 8 x 8 lattices are, 0.6250, 0.7220, and 
0.78125, respectively. Notations for the symbols are summa- 
rized in the legends. 

the wave number of the HOL is (7r,7r/4) and the wave 
number of the LUL is (7r/2, 7r/2). The energy separation 
Ae = e(,r/2,Tr/2) ~ e(ir,w/4) IS a monotonically increasing 
function of t" . We plot the errors % of the ACP approx- 
imation [—Eacp + Eqmc) I Eqmc X 100 as a function of 
Ae in Fig. ^. The error is almost independent of Ae and 
is always less than 1%, as long as the CCM equations con- 
verge. The CCD equations with the ACP approximation 
do not converge when Ae < 0.18. This Ae dependence of 
the errors indicates the enhancements of the errors in the 
larger lattices observed in Fig. |3| comes from the differ- 
ence of the electron fillings p. The CCM is less accurate 
close to the half- filling. 

While it seems the ACP approximation works nicely in 
the Hubbard model up to the intermediate U/t region, 
it fails in the large U/t region and/or close to half- filling 
in the two dimension. The U/t dependence of the error 
was also observed in one-dimensional Hubbard model on 
small clusters. Direct incorporation of r„,n > 3 may be 
necessary in the large U/t region and close to the half- 
filling. The point in this article is that the accuracy of the 
CCM does not depend much on the LUL-HOL gap Ae. 
When we apply the crystal orbital method to solid elec- 
trons, we use a finite number of k points for k space in- 
tegrations and constructing excited configurations. [p|jT^] 
Hence, the Fermi degeneracy in metals is replaced by 
a pseudo-degeneracy among some of the k points. The 
situation is similar to atomic cluster calculations. In nu- 
merical calculations, there are implicitly gaps due to the 
finite number of the k points, even for metallic electrons, 



3 



10- 



o 
o 



+ 



I 



0.0 



0.2 



0.4 
Ae = e 



0.6 



0.8 



1.0 



(II/2.II/2) (lt,ir/4) 



FIG. 4. The error (%) {-Eacp + Eqmc)/Eqmc x 100 
as a function of Ae in the generalized 2D Hubbard model on 
the 8x8 lattice. 



as well as for electrons in semiconductors and insulators. 
The fact that the accuracy of the CCM does not depend 
on such gaps is very favorable, unless numerical results 
are very dependent on details of the calculation param- 
eters such as increments of k rather than on physical 
parameters. Hence, the CCM may be successfully ap- 
plied to ab-initio band calculations on metals as well as 
semiconductors and insulators. 

In conclusion, we have studied the electron correlation 
problem in the doped 2D Hubbard model with the ACP 
approximation to investigate whether or not the approxi- 
mation can be applied to correct the independent particle 
approximations really used in ab-initio band calculations. 
We found most of correlation energies of the 2D Hubbard 
model in the weak and the intermediate U /t region can be 
accounted for by the ACP approximations. In the large 
U /t region and/or close to the half-filling, the error of the 
ACP approximation is not negligible. Further incorpo- 
ration of the electron correlation effect may be necessary 
there. The accuracy of the CCM does not depend on the 
LUL-HOL gap. Hence, the CCM may be successfully ap- 
plied to ab-initio band calculations on metals as well as 
semiconductors and insulators. 
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